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We study the interplay between the Hofstadter butterfly, strong interactions and Zeeman field within the 
mean-field Bogoliubov-de Gennes theory in real space, and explore the ground states of the attractive single¬ 
band Hofstadter-Hubbard Hamiltonian on a square lattice, including the exotic possibility of imbalanced vector 
potentials. We find that the cooperation between the vector potential and superfluid order breaks the spatial 
symmetry of the system, and flourish stripe-ordered Fulde-Ferrell-Larkin-Ovchinnikov (FFLO)-like superfluid 
and supersolid phases that can be distinguished and characterized according to their coexisting pair-density 
(PDW), charge-density (CDW) and spin-density (SDW) wave orders. We also discuss confined systems and 
comment on the likelihood of observing such stripe-ordered phases by loading neutral atomic Fermi gases on 
laser-induced optical lattices under laser-generated artificial gauge fields. 

PACS numbers: 03.75.Ss, 03.75.Hh, 67.85.-Lm, 67.85.-d, 67.80.kb 


I. INTRODUCTION 

The exact energy spectrum of a single quantum particle that 
is confined to move on a two-dimensional tight-binding peri¬ 
odic lattice under the influence of a uniform magnetic flux 
has been known for a long time Qi, where the competi¬ 
tion between the lattice spacing and cyclotron radius gives rise 
to a self-similar complex pattern of sub-bands and mini-gaps. 
However, regardless of all efforts since the prediction of this 
Hofstadter spectrum, there has been very recent but still lim¬ 
ited success in observing some of its signatures in graphene- 
based solid-state materials with artificially-engineered super¬ 
lattices under real magnetic fields M.' In addition, thanks 
to the recent realisation of artificial gauge fields in atomic sys¬ 
tems iMl, there is also an increasing interest on this sub¬ 
ject from the cold-atom community 111^2(11] . In particular, 
by engineering spatially-dependent complex tunneling ampli¬ 
tudes with laser-assisted tunneling and a potential energy gra¬ 
dient, two groups have recently reported realisation of the 
Hofstadter-Harper Hamiltonian using neutral rubidium atoms 
that are loaded into laser-induced periodic potentials |fl6l - [l^ . 

Even though the Hofstadter and Hubbard Hamiltonians 
have themselves been the subject of many works in the 
literature, there has been a lack of interest in the com¬ 
bined Hofstadter-Hubbard Hamiltonian even at the mean- 
field level. For instance, while the use of momentum-space 
BCS formalism limits previous analysis of the attractive 
Hofstadter-Hubbard model only to vortex lattice (VL) config¬ 
urations the existence of pair-density wave (PDW) and 
VL orders have been proposed in the context of a somewhat 
related model; an anisotropic 3D continuum Fermi gas expe¬ 
riencing a uniform magnetic flux ll^ . By first limiting their 
description to the lowest-Landau-level limit and then making 
further assumptions about the strength of the anisotropic trap, 
the authors obtain an effectively a ID Hamiltonian in mo¬ 
mentum space, and solved it using the BCS formalism. The 
existence and characterisation of a variety of distinct stripe- 
ordered many-body phases have either been overlooked or 
gone unnoticed until very recently 12^ . distinguishing our 
work from the literature. 

In particular, here we use Bogoliubov-de Gennes (BdG) 


theory in real space and study the mean-field ground states of 
the attractive single-band Hofstadter-Hubbard Hamiltonian on 
a square lattice, including the effects of imbalanced chemical 
and vector potentials. We find that the cooperation between 
the vector potentials and interaction breaks the spatial sym¬ 
metry of the system, leading to various stripe-ordered super¬ 
fluid (SF) and supersolid (SS) phases that can be distinguished 
and characterized according to their coexisting PDW, charge- 
density (CDW) and spin-density (SDW) wave orders. We also 
discuss possible observation of such stripe-ordered phases by 
confining neutral atomic Fermi gases in laser-induced optical 
lattices under laser-generated artificial gauge fields. 

The rest of this paper is organised as follows. In Sec. [Ill 
first we introduce the physical setting of the problem and 
the model Hamiltonian used, then review the non-interacting 
Hofstadter Hamiltonian and its well-known Hofstadter spec¬ 
trum, and then describe the self-consistent BdG formalism 
which takes fermion-fermion interactions into account within 
the mean-field approximation for pairing. The resultant BdG 
equations are solved in Sec. imi where first we tabulate the 
numerically obtained mean-field ground states, paying a spe¬ 
cial attention to the striped phases in the dimer-BEC limit, 
and then construct the thermodynamic phase diagrams. The 
effects of Hartree shifts on the possible ground states are 
discussed in Sec.|V]in the context of harmonically-confined 
atomic systems. We end the paper with a briery summary of 
our conclusions and an outlook in Sec. [V!] and an Appendix 
comparing the dimer-BEC limit in the Landau and symmetric 
gauges. 


II. THEORETICAL ERAMEWORK 

To explore the ground states of the single-band Hofstadter- 
Hubbard model, we start with 

H = — tijrrola^jrr ~ X! 

ijcr ia 


( 1 ) 
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and consider both thermodynamic and confined systems. 
Here, aj^ {ain) creates (annihilates) a a = {ti-l-} fermion 
on site i, is its hopping parameter from site i to j, and 
= fj. — Vi + h and = fj, — Vi — h are effectively their 
local chemical potentials in the presence of confining poten¬ 
tial Vi and an out-of-plane Zeeman field h. We assume h > 0 
without loosing generality, since h < 0 results can be easily 
deduced by letting The density-density inter¬ 

action term is taken to be local (on-site) and attractive with 
strength g > 0, and the resultant many-body phases are inves¬ 
tigated within the mean-field approximation for the Cooper 
pairs and their superfluidity, as described below. 


largest magnetic field i3-|- = ~ lOOT that is attainable in 

a laboratory, (j>^ and (j>i may be tuned at will in atomic optical 
lattices. 

Let us first set p = 0 and /iio- = 0 in Eq. (|2]), and re¬ 
view the well-known single-particle problem, i.e, the Hofs- 
tadter Hamiltonian for a uniform square lattice. 

B. Hofstadter Butterfly (HB) 

In the non-interacting limit, the single-particle Hofstadter- 
Hamiltonian describing a a fermion can be written as. 


A. Mean-Field Hofstadter-Hubbard Hamiltonian 

In particular, we analyse the following mean-field Hamilto¬ 
nian for square lattices. 


Hofj — / ' I cttj 


CT 7 . ( ^nma-^'n+l.miT 

*2-^‘'"aLaa„,™+i..+H.c. ), (3) 


Hmf — / . 7 , gjqO'ifjO'ia 

ij<7 i<7 

where pi-j- = gi^ — gnu and — gui^ take the Hartree 

shifts into account. Here, riio- = (aj^aia) is the average num¬ 
ber of tr fermions where (• • •) is a thermal average, and the 
remaining terms in Eq. Q involve the complex SE order pa¬ 
rameter Ai = g{ai-faii). These average quantities are spec¬ 
ified below in Eqs. 0-®, and we use them in Sec. nil] for 
characterising the state of the system. 

The complex hopping parameters are assumed to connect 
only the nearest-neighbor sites, i.e, Ljo- = where 

the amplitudes = t > 0 are taken to be equal for 

i and j nearest neighbors and 0 otherwise. The phase, how¬ 
ever, depends on the external magnetic (or artificial gauge) 
field experienced by a fermions. In particular, we use the 
Peierls substitution and take Bij^r = {i/4’0) -^cr(r) • dr, 

with 4>o = ‘l-nh/e the magnetic flux quantum and ACT(r) the 
vector potential which is assumed to be independently con¬ 
trollable for 7 and fermions. Note that while independent 
control of A„ (r) is not possible for conventional solid-state 
materials with real magnetic fields where a corresponds to the 
± projections of spin angular momentum of electrons, such 
a control can be achieved with neutral atomic systems under 
the influence of laser-generated artificial gauge fields where 
pseudo-spin a is just a label for two of the hyperfine states 
of a particular atom. In this paper, we choose Landau gauge 
for the vector potential, i.e, Acr(r) = (0, B^x, 0), leading to 
a uniform magnetic flux = B^i"^ per unit cell penetrat¬ 
ing our square lattice, where ^ is the lattice spacing. Denoting 
{x, y) coordinates of site i by (nf, m£), this gauge simply im¬ 
plies 9ijcr = 0 and 9ij,j = ±2Trn4)a for links along the x and 
y directions, respectively, where 4>a = ^a/ (27r^o) character¬ 
izes the competition between £ and the magnetic length scale 
(cyclotron radius) £b„ — i/fi./ {eB„). We note that while 
</'t ~ '/'I ^ 1 typical electronic crystals, even for the 


where H.c. is the Hermitian conjugate. Eor rational values 
of 4>a = Pal<la, where Per and are positive integers with 
no-common factor, i.e, co-prime numbers, while main¬ 
tains its translational invariance in the y direction, it requires 
qa sites for translational invariance in the x direction. Thanks 
to the Bloch theorem, the 1st magnetic Brillouin zone is de¬ 
termined by —TT < ky£ < TT and —n/q^ < kx£ < 'n'/qa, and 
this increased periodicity motivates us to work with a super¬ 
cell of 1 X sites. The excitation spectrum is determined 
by solving the Schrodinger equation i/oo-T'er = £:{4>(T)'i’<7 for 
all momentum k = {k^, ky) values in the 1st magnetic Bril¬ 
louin zone. Denoting the components of the wave function as 
T'cr = (V'l J V'J J V'a) • ■ •) where 'i/n corresponds to the 

nth site of the supercell, the x Qa Hamiltonian matrix at a 
given k value 


-Cia t; 0 . . 0 r, ■ 

T,y C2a t: 0 . . 0 

0 T, Csa ■ ■ 

. 0 

0 . ... Cn-I.a t; 

T: 0 . . 0 T,, Cna 


(4) 


describes the supercell with periodic Bloch boundary condi¬ 
tions. Here, Cna = —2tcr cos{ky£ + ^-Knya/ga) and = 

The eigenvalues e{4)a) of this g^ x q„ matrix can be nu¬ 
merically obtained for any given rational number 6^ and the 
energy spectrum e{4)a) vs. 4> is known as HB The 

spectrum is shown in Eig.[T| where, for a given 4>a, it consists 
of non-overlapping g^- bands with gcr + 1 energy gaps in be¬ 
tween, and each one of these g^ bands can accommodate 1 /go- 
particle filling with a total filling of 1. Therefore, if we index 
energy gaps as za = {0,1, 2, • • • , go-}, starting from the bot¬ 
tom edge of the band in such a way that the lowest (za = 0) 
and highest (z^ = qa) gaps correspond, respectively, to a par¬ 
ticle vacuum and a fully-filled band insulator, particle fillings 
within all of these gapped regions can be compactly written 
as Za/qcr- Note that while all gaps are open for odd g^, the 
middle Za = gcr/2 gap corresponding to a half-filled lattice 
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FIG. 1. (Color online) The Hofstadter spectrum e/io- is presented as 
a function of cpa- = pa/qa, showing its fractal structure with numer¬ 
ous sub-bands and mini-gaps. 


is not open when is even, and therefore, a half-filled lat¬ 
tice is not an insulator for any q^. In the (p^ ^ 0 limit, the 
HB spectrum recovers the usual tight-binding dispersion of 
cosines = — 2 fo. [cos{kx£) + cos{kyi)] , which has an en¬ 
ergy bandwidth 

Since £ and £b„ are the only two length scales in Eq. (|3ll 
such that (pa = the fractal structure of HB is 

clearly a result of their competition. In addition, each k state 
is ^cr-fold degenerate in the 1st magnetic Brillouin zone (not 
explicitly shown in the figure), i.e, 

^b.kxkyiP^tj) — ^b,kx,ky+27T4>xf (^) 

with b = 1,2, - ■■ ,qa labelling the bands and / = 

1, 2, • ■ • ,qa labelling the degenerate k states. We have re¬ 
cently shown that the HB spectrum plays a crucial role in de¬ 
termining the many-body states of the interacting system 
and our primary objective here is to extend and generalise the 
analysis to imbalanced gauge fields. 


C. Bogoliubov-de Gennes (BdG) Theory 


For this purpose, we diagonalise Eq. (|2]l via 
the Bogoliubov-Valatin transformation, i.e, aia — 
~ where (Tmcr) Creates 

(annihilates) a pseudo-spin cr quasiparticle with energy 
and wave functions Umia and Vmia, and = -Fl and 
S 4 , = —1. The resultant BdG equations can be compactly 
written as. 



£ij'\ 


^ibij 




(7 

mi 1 


where 5ij is the Kronecker delta, and (pl^^ = 

and -UmiiV are the corresponding eigen¬ 

functions for > 0 eigenvalue. Note that the BdG equa¬ 
tions are invariant under the transformation Vmi^ —f Wmit’ 


Umii ->■ ^ therefore, it is suffi¬ 

cient to solve only for Umi = Vmi = and Cm = 
as long as all solutions with positive and negative e™ are kept. 

Using the transformation, the complex order parameter 
can be written as 

Ai = g 'y ( 7 ) 

m 

where/(a:) = l/[e“AfeBT) _|_ jg the Fermi function with ks 
the Boltzmann constant and T the temperature. Equations (| 6 ]l 
and (| 2 ll have to be solved self-consistently for a given p and 
h, such that the total number of a fermions satisfies Na = 
J2i n-icr- Here, 0 < rua = (ol^aia) < 1 is the average number 
of a fermions on site i, and using the transformation, it can be 
written as 


^ ^ l^mzl 

m 

(8) 

'il — ^ ^ /( £771); 

( 9 ) 


m 


for the t and j, fermions, respectively. We note that unlike the 
continuum models where the solutions of the self-consistency 
equations depend explicitly on the high-momentum cut-off, 
requiring a high-energy regularisation in order to obtain cut¬ 
off independent results, the lattice versions given in Eqs. ( 0 - 
(|9]l do not require such a regularisation, since the lattice spac¬ 
ing £ already provides an implicit short-distance cut-off. 

In the absence of gauge fields when = 0, it is generally 
accepted that the mean-field description given above provides 
qualitative understanding either at low temperatures (T <|; Tc) 
for any g or for weak 5 < kU at any T, where Tc is the critical 
SF transition temperature. It is also known that single-band 
Hubbard models gradually become inadequate in describing 
strongly-interacting cold-atom systems on optical lattices, re¬ 
quiring multi-band models ll24l] . In addition, the real-space 
BdG theory goes beyond the standard local-density approxi¬ 
mation since it includes both 0 ^ and Vi exactly into the mean- 
field theory without relying on further approximations. Hop¬ 
ing to shed light on the qualitative effects of gauge fields on 
the ground states of Eq. 0, here we mainly concentrate on 
weak and intermediate 5 at T = 0 as discussed next. 


III. NUMERICAL FRAMEWORK 


In order to explore the possible phases, let us set U = 0 and 
consider a uniform 45f x 45f square lattice, which is large 
enough to construct the thermodynamic phase diagrams for 
pa — {0,1/6, ±1/4}. We note that even though our phase 
diagrams are reliable, the phase boundaries should be taken 
as qualitative guides to the eye due to the possibility of minor 
finite-size effects. We neglect the Hartree shifts for the mo¬ 
ment because not only the self-consistent solutions converge 
much faster but also the resultant phase diagrams are much 
more easier to interpret and understand. In addition, since 
none of the PDW, CDW and SDW instabilities are driven by 
these shifts, our qualitative mean-field results already paves 
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TABLE I. While the S-SS* phase has a small hut finite sign-changing 
striped-SDW order, the system is globally unpolarized very much 
like the unpolarized uniform superfluid (U-SF) or unpolarized striped 
supersolid (S-SS) phase. 


Phase 

|A.| 

Tii^ + 


0(T 

U-SF 

Uniform 

Uniform 

0 

Pr = —pi 

S-SF 

PDW 

1 

0 

II 

•<- 

S-SS 

PDW 

CDW 

0 

Pi = Pi 

S-SS* 

PDW 

CDW 

SDW 

\pi\ 7^ \pi\ 

P-SF 


otherwise 



the way to quantitative understanding of the possible ground 
states of Eq. 0 . However, see Sec. IVBI for the effects of 
Hartree shifts on confined systems. 

For this purpose, we numerically solve Eqs. ©-dill at T = 
0 , and obtain self-consistent solutions of A^/f and riia as 
functions of g/t, g/t, h/t and (pa-. This can be achieved nu¬ 
merically via the iterative method of relaxation as follows. For 
a given set of parameters, first (i) start with an input set of A^, 
then (ii) construct the BdG matrix given in Eq. (|6]l, and then 
(iii) use its eigenstates in Eq. 0 to generate a new set of A^, 
and finally (iv) repeat these steps until the input and output 
sets of Ai lie within a confidence level. Once this iterative 
method converges, (v) use Eqs. 0-0 to calculate riia- It 
turns out that while Eqs. @-0 have unique solutions in the 
low-/i/p limit, they in general allow for multiple solutions for 
the polarized many-body phases, and therefore, it is essential 
to try several initial sets of A^ and verify the (meta)stability 
of the solutions. 


A. Ground-State Phases 

Depending on the spatial profiles of |Ai|, rii^ and riii, 
we distinguish the single-particle band insulator and normal 
phases from the ordered many-body ones using the following 
criteria. When h/g is sufficiently high that A^ —0 (pre¬ 
cisely speaking |A*| < 10 in our numerics) for every i, 
the ground state can be a cr-vac phase which is a vacuum of 
a component with riia = 0 , a a-I{m/n) phase which is a 
band insulator of a component with uniform Uia = m/n, 
a a-N phase which is a normal cr component, or an p'P-PN 
phase which is a polarized normal mixture of t and I com¬ 
ponents. We checked in our numerics that while a-N and 
pl-PN phases have slightly non-uniform riia for (p ^ 0, 
the C 4 symmetry of the square lattice is preserved. On the 
other hand, when h/g is sufficiently low that A^ 7 ^ 0 (i.e, 
|A*| > 10 ^t) for some i, the ground states can be character¬ 
ized according to Table U Unlike our earlier work ll2^ . here 
we do not finely classify the polarized superfluid (P-SF) phase 
depending on the coexisting (striped or non-striped) PDW, 
CDW, SDW and/or VL orders. Instead, we focus mostly on 
the existence of striped phases in the dimer-BEC limit as the 
main message of this manuscript, for which physical (analyt¬ 
ical) insight are also given. 

The globally-unpolarized states are denoted by U-SF, S- 



FIG. 2. (Color online) Characterisation of globally-unpolarized 
many-body phases, (a) Typical |Ai|/i profiles are shown for the 
U-SF (left) and S-SF phases (right), where pf = = 1/4 and 

p-f = Pi = i- /4, respectively, and p = 0 (uniformly half-filled) in 
both figures, (b) Typical |Ai|/f (left) and rii-f -f (right) profiles 
are shown for the S-SS phase, where pi = pi = 1/4 and g = —t. 
(c) Typical |Ai|/f and rii-f — Uii profiles are shown for the S-SS* 
phase, where p.f = 0, pi = 1/4 and g — —t. Note in (c) that 
even though the system is globally unpolarized, it has both SDW and 
CDW (not shown) orders. Here, {x,y) are in units of £, and we set 
h — 0 and p = 7f in all figures. 


SF and S-SS, and they stand, respectively, for uniform-SF, 
striped-SF, and striped-SS. The S-SS* state is also globally 
unpolarized, very much like the S-SS phase but it has an ad¬ 
ditional sign-changing striped-SDW order driven by the im¬ 
balance between and pi. For instance, typical |Ai| and 
riii ± riii profiles are illustrated in Fig.|2]for all of them. De- 
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pending on p,, h, and one of the U-SF, S-SF, S-SS and 
S-SS* phases always appears in the thermodynamic phase di¬ 
agrams beyond a critical g/t threshold, as discussed next. 

B. Dimer-BEC Limit in the Landau Gauge 

When g/t ^ 1 is sufficiently high, the physics must 
eventually be determined by the two-body bound states, i.e. 
Cooper pairs become bosonic dimers, and unless g/t ^ oo, 
the dimer-dimer interaction [gdd ^ is finite. 

Such weakly-repulsive dimers can effectively be described 
by the Hofstadter-Bose-Hubbard model, where superfluidity 
has recently been shown to break translation symmetry in the 
weakly-interacting limit ll^ . 

In the ideal-dimer limit of our model Hamiltonian, the 
only way a tightly-bound dimer to move from a site i to j 
in the lattice is via what is known as pair-breaking mecha¬ 
nism, i.e, virtual ionisation of its constituents costs a penalty 
of g, and this gives rise to the effective dimer hopping pa¬ 
rameter tijd = 2tij^tij^/g. Therefore, the effective hop¬ 
ping amplitude and gauge field of the dimers can be written 
as td ~ ^t^ti/g and pd = Pt + h = Pd/<ld, respectively, 
where pd = {p^q^ +Viqt)/Q ^d = q^q\./Q- Here, Q is 
a positive integer number chosen such that pd and qd are co¬ 
prime numbers, and it depends on the entire {pt,P4., (Jti 94.} 
set. Since HB for dimers is g^-fold degenerate, their BEC 
order parameter has contributions from all degenerate kd = 
{( 0 , 0 ); (0,27r(/)d//f)} momenta, where / = 1 , • ■ ■ , - 1 

such that T'id = cq -f and c/ = |c/|e*'’-t 

are complex variational parameters. However, unlike atomic 
bosons where all of the degenerate states have equal weight, 
dimer bosons are fermion pairs and the number of ways of 
creating them with kyd = ky^ -\- kyi momentum depends on 
/, (j>^ and (pi. For instance, there are 2{q — f) — 1 ways of 
intra-band pairing when = 4>i = p/q and q is even. Thus, 
this analysis show that higher kyd states contribute less and 
less, forming a perturbative series. 

It turns out that the first order (/ = 1 ) correction is already 
much smaller than the zeroth order (/ = 0 ) one, and that the 
/ > 2 terms are always negligible when g/t is sufficiently 
large. This is because all of our numerical results fit quite 
well with 

|Ai| = |Ao| -f |Ai| [1 - cos{2^T(()div/^)] , (10) 

in the entire globally-unpolarized region, including S-SF, S- 
SS and S-SS* phases. Here, the kd = (0,0) contribution 
|Ao| = {g/2 — /g')^n(/l — n) is uniform in space and 

determined by the total average filling n with p = (p /2 — 

— 1) i^ . |Ai| « t^/g for p Ri 0, and iy is 
the y coordinate of site i. Moving towards the BCS side, 
the second-order correction to Eq. (fTOl i can be shown to be 
-I-IA 2 I cos{A'K(j>diy/(-) for even qd- Since this term is in- (out- 
of-TT-) phase with the zeroth (first) order term, it tends to open 
throughs along the peaks arised from the first order one, sug¬ 
gesting that ttf — tto = tt/, i.e., the form of coincides 
with I Ai I under these conditions. Equation ( fTOl l clearly shows 
that modulations of | A^ | have a spatial period of qd lattice sites 


along the y direction. It also implies that it is the cooperation 
between pd and g that is responsible for the broken spatial 
symmetry and appearance of stripe order, and even though the 
stripe order gradually fades away with increasing g, it survives 
even in the g ^ W limit as long as g/t is finite. 

Thus, this analysis suggests that the existence of stripe- 
ordered SE and SS phases is not an artefact of the mean-field 
description, and they are physically expected in the dimer- 
BEC limit of the attractive Hofstadter-Hubbard model, as dis¬ 
cussed next. 


IV. THERMODYNAMIC PHASE DIAGRAMS 

Despite tremendous efforts over several decades, while the 
exact phase diagram of even the simplest Hubbard model 
(which does not include the gauge fields or Zeeman fields) 
is still the subject of a hot debate, the mean-field phases and 
resultant phase diagrams of the mean-field Hubbard model are 
pretty much settled. To appreciate the effects of gauge fields, 
first we study Eq. © with = 0 . 

A. No Gauge Fields: p^ = pi = 0 

Our results for this limit is presented in Eig. [3] where we 
set p, = 0 in [ 2 ^ a) corresponding to a half-filled lattice, and 
p = —t in[3b). We find that the phase diagrams are very sim¬ 
ilar, and depending on the particular value of g, there are two 
critical h fields. Since EELO phase occupies a tiny parameter 
space near the normal phase boundary and only on the BCS 
side when g/t < W, we do not finely classify the character 
of P-SE phase in Eig. [3 and throughout this paper. The U-SE 
phase, where A^ = Aq for all i, turns into a P-SE beyond a 
first critical field hc^ , and then the P-SE phase becomes an 
PN beyond a second critical field hc2 > hci . Our numerical 
results indicate that hc^ ^ |Ao| where |Ao| is evaluated at 
h = 0 for the same parameters. 

In the strongly-interacting limit when p S> f, it can be 
analytically shown for thermodynamic systems that |Ao| = 
{g/2 — At^/g)^n{2 — n), where n = n-i- -f nj, is the total 
fermion filling. We checked that this thermodynamic expres¬ 
sion agrees very well with our finite-lattice results, as it gives 
|Ao| ~ 7.23f for p = 0 or n = 1 and |Ao| ~ 7.18f for 
p = —f or n Ri 0.875 when g = 15f, while we find, respec¬ 
tively, |Ao| Ri 7.25f and |Ao| ~ 7.19f for the same parame¬ 
ters in our BdG calculations. In the weakly-interacting limit 
when g is sufficiently small so that Ai 0 for every i, we 
note that the system will be a {.-vac for h > At when p = 0 
and for h > 3t when p = —t. Next, we are ready to discuss 
the effects of balanced gauge fields. 

B. Balanced Gauge Fields: pf- = pi ^ 0 

In Eig. m we present the p^r = 1/4 phase diagrams for 
p = 0 inUa) and p = —t in|4}b). The p — 0 case is very spe¬ 
cial since it corresponds to a half-filled lattice with particle- 
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FIG. 3. (Color online) No gauge field: cf)^ = (f)^ = Q case. The 
ground-state phase diagrams are shown for /r = 0 in (a) and fi = —t 
in (h). 


hole symmetry, where rii-f + riii — 1 independently of i, no 
matter what the rest of the parameters are. In comparison to 
Fig. [3] the t/fg. = 1/4 diagrams have much richer structure 
involving large regions of stripe-ordered phases. To under¬ 
stand the physical origin of the resultant phase diagrams and 
stripe order, next we discuss the analytically tractable high- 
and low-Zi/g limits. 

When h/g is sufficiently high, we can directly read off the 
single-particle ground state of the a component from HB for 
any given For (pa- = 1/4, the energy spectrum consists of 
4 bands: the tr component is a a-vac for /Xcr < — 2.83f, a a-N 
for -2.83f < fia < -2.61f, a cr-/(l/4) for -2.61t < Ha < 
-1.082t, a a-N for -1.082f < < 1.082f, a cr-/(3/4) for 

1.082f g-a PP, 2.61t, a a-N for 2.61f PP fJ-a Pp, 2.83f and a 
tT-/(l/l) for 2.83t < /icr. Using = g,-\- h and /rj, = /i — ft. 
in these expressions, the high-ft/p structure of Fig. |4] imme¬ 
diately follows. As h/g gets smaller, the single-particle I 
and N phases must pave the way to ordered many-body ones, 
as increasing the strength of the pairing (attractive potential) 
energy eventually makes them energetically less favourable. 
For pa = 0, it is intuitively expected and numerically con¬ 
firmed above that the fi-UTV to P-SF phase transition bound¬ 
ary p(ftc) is a monotonic function of ft, which is simply be¬ 
cause the non-interacting system has a very simple band struc¬ 
ture with cosine dispersions. However, due to the presence 
of multiple bands, the transition boundary (?(ftc) becomes a 
complicated function of ft for finite pa- For instance, we find 
a sizeable hump in Fig.Ua) around ft k, 2.It and another one 
in Fig.Ub) around ft ss 1.7t, the peak locations of which co¬ 
incide intuitively with the pP-PN regions that are sandwiched 
between VAC and/or I. 

On the other hand, when h/g is sufficiently small, the 
ground state is expected to be an ordered many-body phase 
with no polarisation. In sharp contrast to the pa — 0 case 




g/t 

FIG. 4. (Color online) p-^ = p_i — 1/4 case. The ground-state phase 
diagrams are shown for p = 0 in (a) and g = —tin (b). 


where U-SF is numerically confirmed above to be the ground 
state for any g, we show in Fig. |4]that S-SF and S-SS are, 
respectively, stable for g = 0 and g = —t when pa = 1/4. 
Note that since g = 0 corresponds to half filling for any pa, 
the unpolarized ground states necessarily have uniform fill¬ 
ings, i.e, rii^ = riii = 1/2 for every i. Therefore, in the 
low-ft/p limit, while only | | is allowed to have spatial mod¬ 
ulations in Fig.UJa), both | | and riia modulates in Fig.|4jb). 

In comparison, the pa = 1/6 phase diagrams are shown 
in Fig. |5] and they are in many ways similar to the pa = 1/4 
ones. The main difference is in the high-ft/p limit which again 
directly follows from HB. For pa = 1/6, the energy spectrum 
consists of 6 bands: the a component is a cr-vac for ga < 
—3.076f, a a-N for a narrow band around ga ~ —3.076/, a 
cr-/(l/ 6 ) for -3.076/ < ga < -1.59/, a a-N for -1.59/ < 
ga PP, —1.41/, a cr-/(l/3) for —1.41/ ^ ga PP —0.65/, a a-N 
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g/t 


FIG. 5. (Color online) (j>^ = = 1/Q case. The ground-state phase 

diagrams are shown for ^ = 0 in (a) and /i = —f in (b). Note in (b) 
the presence of a sliver of \.-N region around h « 2.076f. 


for —0.65f < ficr < 0.65f, a ct-/( 2/3) for 0.65f 
lAlt, a a-N for 1.41f < /io- < 1.59f, a cr-/(5/6) for 1.59f < 
Mo- ^ 3.076f, a a-N for a narrow band around ?a 3.076f 
and a-I{l/l) for 3.076f < ^o-. As a consequence of this, we 
note in Fig. [S^b) that the system intuitively requires a hnite 
threshold for g/t even at h = 0 , in order to develop any kind 
of many-body order. In addition, it is intriguing to see that 
the sliver of \.-N region that is sandwiched between 4,-vac and 
4,-/(l/6) around h Ri 2.076f gives rise to a sizeable hump in 
Fig. I^b). This is clearly a result of increased single-particle 
density of states. 

Note in Figs. [35] that the transition from an unpolarized to 
a polarized ordered phase occurs at a lower h for any given g 
as is increased from 0. This is a consequence of smaller 
non-interacting energy bandwidths: as (/>„ increases from 0 to 


1/6 to 1/4 then W shrinks from 8 t to 6.15f to 5.65f, mak¬ 
ing it possible to polarize the ground state with a smaller 
and smaller h. In Figs. |4| and 0 the P-SF regions are dom¬ 
inated mainly by a phase that can be characterized by almost- 
striped PDW and SDW orders with some additional corruga¬ 
tions along the stripes that is caused by h 7 ^ 0. For instance, 
when this phase is nearby to an insulating one, it generally has 
a very small SDW order in the background on top of a large 
and uniform polarisation. 


C. Imbalanced Gauge Fields: 7^ ^4. 

As we argued in Secs.Hland lll A1 while independent control 
of the gauge fields (pf and pi is not possible for conventional 
solid-state materials with real magnetic fields, such a control 
is plausible with neutral atomic systems. Motivated by this 
exotic possibility, here we study two different limits. 



FIG. 6. (Color online) Time-reversal symmetric gauge fields: = 

—(f>l = 1/4 case. The ground-state phase diagrams are shown for 
p = 0 in (a) and g = —tin (b). 

As the first limit, we consider a pair of time-reversal sym¬ 
metric gauge fields, i.e, = —pi- For instance, cp^ = 1/4 
phase diagrams are shown in Fig. jbj where we set p = 0 
in|6ja) and g = —t in|6jb). Thanks to the time-reversal sym¬ 
metry, even though the ground state is not a P-SF but an unpo¬ 
larized SF at h = 0, it is not properly indicated in these figures 
for low g/t. The general structures of the transition bound¬ 
aries that are seen in these phase diagrams are quite similar to 
the ones shown in Fig.|4|for the = 1/A case. How¬ 

ever, there is an important caveat in the dimer-BEC limit: the 
ground state becomes a U-SF for any g as long as h/g is suffi¬ 
ciently low. Given our analysis in Sec. lIIIBl this is intuitively 
expected since the effective gauge field of Cooper pairs vanish 
(pd = 0) in the dimer-BEC limit as the gauge field of p and 
4 , fermions precisely cancel each other. In addition, the P-SE 
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regions necessarily shrink here, since the U-SF to P-SF transi¬ 
tion boundaries are expected to be close to the no-gauge-field 
{(j)a = 0) ones shown in Fig. [3 



g/t 

FIG. 7. (Color online) Charged-uncharged mixtures of fermions: 
01- = 0 and 04. = 1/4 case in (a-b) and 0i- = 1/4 and 04, = 0 
case in (c). The ground-state phase diagrams are shown for p = 0 
in (a) and p = —f in (b-c). Thanks to the particle-hole symmetry, 
the 01- = 1/4 and 04. = 0 phase diagram for /r = 0 can easily be 
deduced from (a) via t—and 

As the second limit, we set one of the gauge fields to zero, 
corresponding effectively to a charged-uncharged mixture of 
two-component fermions. For instance, (0-|- = 0, 04 . = 1/4) 
phase diagrams are shown in Fig. [T] where we set /r = 0 in 


iTja) and ^ = —t in|3b), and (0i- = 1 / 4 , 04 . = 0) diagram 
is shown in Fig. [TJc) where we set p = —t. Thanks to the 
particle-hole symmetry around half-hlling, (0i- = 1/4,04^ = 
0 ) phase diagram for ^ = 0 can easily be deduced from|3a) 
via and 0— and therefore, it is not shown. Since this 
symmetry also prevents polarisation at h = 0, even though the 
ground state is not a P-SF but an unpolarized non-uniform (but 
non-striped) SF for weak g/t, this is not properly indicated 
in Fig. 13 a). However, the imbalance between gauge fields 
causes P-SF in Figs.|7Jb) and|3c) even at h = 0. Similar to 
the analysis given in Sec. IIV Bl the high- and \ow-h/g limits 
can be directly read off from HB and effective dimer-BEC de¬ 
scriptions, respectively, with again an important caveat in the 
dimer-BEC limit: the ground state becomes a S-SS* for g ^ 0 
as long as h/g is sufficiently low. As shown in Eig. |2jc), in 
addition to the coexisting striped-PDW and -CDW orders, S- 
SS*has an additional sign-changing striped-SDW order driven 
solely by 0i- ^ 04.. Note also that if (0i- 0, 04 . = 0) then 

all of the coexisting orders of S-SS* phase are periodic along 
the y direction with periodicity qd = since 0^ = 0i-. 


D. Stripe Order vs. FFLO Modulations 

It is clearly the cooperation between g, 0i- and 04^ that is 
responsible for the broken spatial symmetry and appearance 
of stripe order, causing much more prominent stripes for in¬ 
termediate 5 at a given h. The stripe order is a direct result 
of HB: for a given </>„, the spectrum consists of gc-bands in 
the 1st magnetic Brillouin zone within which each k state is 
qcr-fold degenerate. Therefore, when p 0, not only intra- 
and inter-band pairings but also pairings with both 0 and a 
set of non-zero center-of-mass momenta are allowed j^l^ . 
leading to a non-uniform | A,:| with spatially-periodic modu¬ 
lations, e.g, a PDW order 11:^ . The directions of center-of- 
mass momenta determine the direction of modulations, mak¬ 
ing it gauge dependent, e.g, y direction in Fig. |3 When the 
striped-PDW order is sufficiently large, it drives an additional 
striped-CDW order in the total fermion filling, giving rise to 
striped-SS phases. 

We emphasise that the instabilities towards stripe-ordered 
phases discussed in this paper are driven by the gauge fields, 
and they may formally not be identified with the FFLO 
phase which is driven by the Zeeman held and is character¬ 
ized by cosine-like sign-changing I A,; I oscillations along a 
spontaneously-chosen direction 112^4^ . In addition, while 
the periods of our striped-PDW, -CDW and -SDW orders are 
always given by qd, the period of FFLO modulations is de¬ 
termined by the mismatch h between 0 and 0 Fermi surfaces. 
For instance, when 0i- = 04^ = p/q, the stripes have a spa¬ 
tial period of q or q/2 lattice sites, depending on whether q is 
odd or even. Lastly, while our striped phases survive even in 
the extreme dimer-BEC limit (g/t ^ 1) for a large parame¬ 
ter space, the FFLO modulations survive not only in the BCS 
limit but also for a tiny parameter space nearby the P-SF to N 
transition boundary. 
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V. CONFINED ATOMIC SYSTEMS 

Having explored the ground states and phase diagrams of 
thermodynamic systems, here we study confined systems and 
comment on the likelihood of observing stripe-ordered phases 
by loading neutral atomic Fermi gases on laser-induced opti¬ 
cal lattices under laser-generated artificial gauge fields. For 
this purpose, we consider a harmonically-confined 51£ x hli 
square lattice with an isotropic trapping potential Vi = a|rip 
centered at the origin, where a = 0.01t/£^ is its strength and 
Ti = (ix, iy) is the position of site i. 


A. Effects of Harmonic Confinement 

The local ground states of trapped systems can be reliably 
inferred through the so-called local-density approximation, 
where the local density of the system is mapped to that of 
a thermodynamic one with the same density. This description 
is known to be very accurate for large systems that are trapped 
in slowly-varying potentials. For our model Hamiltonian, due 
to the energy gaps of HB and the Pauli exclusion principle, 
one expects the so-called wedding-cake structures in Ui^ and 
Ui^ profiles of non-interacting fermions at T = 0, where the 
number of mini-gaps determines the number of spatially-flat 
riia regions for a given Thus, wedding-cake structures 
consist of a number of insulating regions that are sandwiched 
between normal regions. However, since the majority of these 
mini-gaps are very small compared to t, finite g and/or finite 
T quickly smear out the flat regions, making their detection 
nearly impossible. In sharp contrast, here we show that the 
broken spatial symmetry and stripe orders persist at interme¬ 
diate and strong interactions, providing a viable knob for the 
experimental probe of the fractal structure of HB. 



FIG. 8. (Color online) The trap profiles are shown for (j>^ — 1/4, 
(j>l — 0, g = t, h — 0 and g = 5t. Here, (a;, y) are in units of 1. 

In Fig. [ 8 ] we illustrate a typical self-consistent solution for 
a trapped system when = l/i, fpi = 0, IJ- = t, h — 0 


and g = 5t. The total numbers of a fermions are approxi¬ 
mately given by N-^ = N_i « 464. While the remnants of 
the so-called wedding-cake structure, i.e., spatially-flat rii'^ 
regions around integer multiples of 1/4 fillings, are hardly 
recognisable, a large PDW order is clearly visible. Given 
the phase diagrams discussed in Sec. IIV Cl both CDW and 
SDW orders are expected to be weak around half-filling, since 
ftit + ’^*4 ~ 1 *^he center of the trap for this particular set 
of data. 



y 

FIG. 9. (Color online) The trap profiles are shown for a; = 0 cuts 
along the y direction which is in units of 1. Here, />-[■ = c/j, = 1/4, 
g — t and h = 0, and therefore, the system is locally unpolarized at 
every i. 

It is easier to visualise and present such trap profiles for a 
cut along the y direction at a particular x value. For instance, 
we show a; = 0 cuts in Figs.l9land[T0l where = 1 /A, g = t 
and = 0 in both figures, but = 1/4 and p^ = 0, respec¬ 
tively. While the local ground states are always unpolarized 
in Fig. |9] where ni\ = for every i, the imbalance between 
p^ and Pi causes small but visible SDW orders in Fig.[T0l We 
note that p^ ^ p^ may also cause a global polarisation, i.e, 
Ni 7 ^ Ni, for weak g, however, this polarisation must gradu¬ 
ally disappear towards the dimer-BEC limit. For instance, as 
g/t increases to (4,5, 6, 7), while is approximately 

given by (454,468,491, 519) in Fig.|9l and are given, 
respectively, by (456,464,489, 518) and (448,464,489, 518) 
in Fig.fTOl 

These figures show that the CDW and SDW orders tend to 
be more prominent for intermediate g as long as the system 
is away from half-filling. This is quite intuitive since the ap¬ 
pearance of a PDW order breaks the spatial symmetry of the 
system at the first place. The spatial periods are, respectively, 
given by 2 and 4 sites in Figs. l9]andfT0l and these findings 
are in agreement with our analysis given in Sec. lIIIBl In ad¬ 
dition, since the relative stripes eventually fade away towards 
the dimer-BEC limit, the trap profiles slowly recover the usual 
(no-gauge-field) results in both figures. It is also pleasing to 
see that the valleys of the PDW and CDW orders and peaks 
of the SDW order coincide when they coexist. These results 
suggest that observation of PDW, CDW and SDW features as 
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FIG. 10. (Color online) The trap profiles are shown for a; = 0 cuts 
along the y direction which is in units of i. Here, — 1/4, (j)_i = 0, 
y, = t and h — 0. 

a function of magnetic flux may furnish clearest and direct ev¬ 
idence for the existence of multiple band structure, and hence 
indirectly for the fractal HB, in trapped atomic systems. 


B. Effects of Hartree Shifts 

Since most of our phases have either coexisting CDW 
and/or SDW orders, our phase diagrams may not be conve¬ 
nient to generate more accurate phase diagrams by including 
the Hartree terms via a simple shift in However, we 

still neglected these shifts in our diagrams for their numeri¬ 
cal as well as analytical simplicity. For instance, including 
these shifts in the self-consistency Eqs. (IDl-® not only re¬ 
quires about an order of magnitude more iterations to con¬ 
verge, but also it complicates our current intuition making 
it more difficult to extract the relation between HB and the 
non-monotonic dependences of some of the phase boundaries. 
Note that since Hartree shifts have no role in driving the stripe- 
ordered phases, which is particularly clear in the dimer-BEC 
limit where do not explicitly play any role in our analy¬ 
sis, their inclusion is expected to change some of the transi¬ 
tion boundaries without much effect on the stability of phases. 
Eurthermore, since the mean-field theory provides only a qual¬ 
itative description of the phase diagrams and the accuracy 
of our results can be somewhat improved by including these 
shifts, one still needs to go beyond this approximation for ex¬ 
perimentally more relevant diagrams. Therefore, even though 
Hartree shifts are neglected in Sec. lIVl our results may already 
pave the way to qualitative understanding of the exact ground 
states of the attractive Hofstadter-Hubbard model. 

To illustrate these points, the Hartree-shifted trap pro¬ 
files are shown in Eigs. [TT] and [12] for the parameters of 
Eigs. |9| anddOj respectively. Comparing these figures show 
that while the inclusion of the Hartree shifts does not have 
much effect on |Ai| for these particular sets of data (thanks 
to the particle-hole symmetry around half-filling), it affects 



FIG. 11. (Color online) The trap profiles are shown for x = 0 cuts 
along the y direction which is in units of f. Here, (j>^ = = 1/4, 

y = t and h = 11, i.e., same as Fig.j^with the Hartree shifts included. 


the total filling quite a bit. Eor instance, as g/t increases 
to (4, 5,6, 7), while is approximately given by 

(288,279,279,286) in Eig. [TT] and are given, re¬ 
spectively, by (298,278,278,285) and (276,282,278,285) 
in Eig. [T2| However, the visibility of the striped-PDW and 
-CDW orders remain largely the same in both cases. In ad¬ 
dition, we note that the remnants of the wedding-cake struc¬ 
tures, i.e., spatially-flat ny -F regions around 1/2 fillings, 
are almost recognisable in Eigs. [TTT a) and [T2] a) when g = 4t 
or less (not shown). While the non-interacting / fermions are 
insulating at 1/4 filling in both figures, the non-interacting / 
fermions are insulating (normal) in Eig. [TTIdT^ . Thus, these 
insulating regions leave their traces as distinct |Ai| dips in 
both figures near y = 13£ when g is sufficiently weak. 



FIG. 12. (Color online) The trap profiles are shown for x = 0 cuts 
along the y direction which is in units of £. Here, t/f = 1/4, t/i = 
0, y = t and h = 0, i.e., same as Fig. [T^ with the Hartree shifts 
included. 

Given these numerical illustrations, it is clear that our phase 
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diagrams already shed some light on a new stripe mechanism 
in the dimer-BEC limit, showing that the fate of stripe-ordered 
SF and SS phases are not affected by the Hartree terms. Hav¬ 
ing discussed the effects of conhnement potentials, we are 
ready to end the paper with a briery summary of our conclu¬ 
sions and an outlook. 


VI. CONCLUSIONS 

Our mean-held results for the attractive single-band 
Hofstadter-Hubbard model on a square lattice are as follows. 
In the presence of a Zeeman held h, in addition to the in¬ 
triguing phase transition boundaries between the N/WAC and 
SF phase, we found a number of distinct many-body phases 
which can be characterized with respect to their coexisting 
striped-PDW, -CDW and -SOW orders. Even at h = 0, 
we reached four important conclusions. First, we numeri¬ 
cally found an unpolarized striped-superhuid phase (S-SF) in 
a large parameter space. Unlike the conventional FFLO phase 
which is driven by h, our S-SF is driven only by the gauge 
helds. Second, we numerically found an unpolarized striped- 
supersolid phase (S-SS) in a large parameter space. Unlike 
the conventional SS phase which is yet to be observed and 
is driven either by long-range (e.g., nearest-neighbor) interac¬ 
tions or the presence of a second species (e.g., Bose-Fermi 
or Bose-Bose mixtures), our S-SS is again driven only by 
the gauge helds. Third, we also found a locally polarized 
but globally unpolarized striped-SS phase (S-SS*) when the 
gauge helds are imbalanced. Lastly, we provided analyti¬ 
cal insights on the microscopic origins of these stripe-ordered 
phases, suggesting a new physical mechanism that gives rise 
to FFLO-like SF and SS phases in the dimer BEC limit. 

The importance of these results can be highlighted as fol¬ 
lows. First, spatially-modulated SF and SS phases are both 
of high interest not only to the atomic physics community but 
also to the condensed-matter, nuclear and elementary-particle 
physics communities. Second, the unusual appearance of the 
stripe order is very exotic and fundamentally important by it¬ 
self, because the connection between the striped-charge order 
that is observed in copper-oxide materials and the formation 
of high-Tc superconductivity has been the subject of a long 
debate in the literature. Even though our work offers no di¬ 
rect relation to cuprate superconductors, understanding stripe- 
ordered phases in the cold-atom context may still prove to be 
benehcial for the high-Tc community. Third, the existence of 
stripe-ordered phases is not an artefact of our mean-held BdG 
description, since they are analytically motivated in the dimer- 
BEC limit. Therefore, we highly encourage further research 
in this direction with different lattice geometries, gauge helds, 
etc., in particular the beyond mean-held ones. 
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Appendix A: Dimer-BEC Limit in the Symmetric Gauge 

It may be important to remark here that the analysis given 
in Sec. IIIIBI depends on the particular artihcial gauge held 
that is simulated in a cold-atom experiment. Next, we use a 
rotationally-invariant symmetric gauge for the vector poten¬ 
tial, i.e, Ao-(r) = Bcr{—y,x,0)/2, and analyse the spatial 
structure of the order parameter |Ai| in the dimer-BEC limit, 
which clearly reveals this dependence. 



FIG. 13. (Color online) Typical |Ai|/t prohles are shown for the 
checkerboard-like stripes, where (j>f = (pi = 1/4 in the left and 
1/8 in the right hgure as determined by the symmetric gauge. Here, 
{x, y) are in units of I, and we set p = 0, /i = 0 and g = 8t in both 
hgures. 

Typical |Ai| prohles are illustrated in Fig.[T3]for (pa- = 1/4 
and 1/8, i.e. cpd = + (pi = Pd/qd is 1/2 and 1/4, re¬ 

spectively. Since the C 4 symmetry of the square lattice is pre¬ 
served in this gauge, the stripes are checkerboard-like in the 
(x, y) plane, such that generalisation of Eq. (fTOl) to 

I Ail = |Ao| + I All [2 - cos{'K(pdix/tj - cosiiicpdiy/i)] 

(Al) 

hts very well with all of our numerical results in the dimer- 
BEC limit. Here, the uniform contribution |Ao| = (p/2 — 
jg)\Jn{2 — n) is determined by the total average hlling 
n with y = (p/2 — 8t‘^/g){n — 1) ll^ . |Ai| Ri f^/p 
for /i « 0, and ri = {ix,iy) is the position of site i. We 
again note that this analytical expression is consistent with 
our expectation that while the dimer-BEC order parameter 
in principle has contributions from all degenerate momenta 
4'id = Skdwhere Ckd = |ckd |e*'’‘‘'i are com¬ 
plex variational parameters, only the lowest order kd = 
{{0,0)','^(pd{fx,fy)/i} terms contribute when g/t is suffi¬ 
ciently large. This is because dimer bosons are fermion pairs 
and the number of ways of creating them with kd = -F k^ 
momentum depends on fx, fy, (p^ and (p^^, and higher kd 
states contribute less and less, forming again a perturbative 
series. 

In addition, moving towards the BCS side, the 
second-order corrections to Eq. (lAlT l can be shown 
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to be -\-\I^ 2 \[cos{ 2 'K(j)dix/(-) + cos(27r^d*y/-^)] and 
+ IA 2 I cos(7r(/)dia;/^) cos(7r^d*i//^) for any qd- Since 
these terms are in- (out-of-7r-) phase with the zeroth (first) or¬ 
der term, they tend to create dimples at the peaks arised from 
the first order one, suggesting that — 1^0 = 7r(/a; -1- fy), 


i.e., the form of again coincides with |Ai| under 

these conditions. Finally, we note that both our numerical 
results as well as the analytical fit Eq. (lAlt clearly show that 
modulations of |Ai| have a spatial period of 2 qd lattice sites 
along both x and y directions, as expected for the resultant 
dimers in the symmetric gauge. 
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